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We study the non-equilibrium slow dynamics for the Kitaev model both in the presence and the 
absence of disorder. For the case without disorder, we demonstrate, via an exact solution, that 
the model provides an example of a system with an anisotropic critical point and exhibits unusual 
scaling of defect density n and residual energy Q for a slow linear quench. We provide a general 
expression for the scaling of n (Q) generated during a slow power-law dynamics, characterized by 
a rate r _1 and exponent a, from a gapped phase to an anisotropic quantum critical point in d 
dimensions, for which the energy gap Ag ~ fcf for m momentum components (i = l..m) and 
~ kf for the rest d-m components (i = m + l..d) with z < z : n ~ T -l™.+(d-m.)z/z']ua/(zv a +i) 
(Q ~ r -[(m+*)+(d-«0*/*V«/(*''«+i)). Thege general expressions reproduce both the corresponding 
results for the Kitaev model as a special case for d — z' = 2 and m — z = v — 1 and the 
well-known scaling laws of n and Q for isotropic critical points for z — z' . We also present an 
exact computation of all non-zero, independent, multispin correlation functions of the Kitaev model 
for such a quench and discuss their spatial dependence. For the disordered Kitaev model, where 
the disorder is introduced via random choice of the link variables D n in the model's Fermionic 
representation, we find that n ~ r _1//2 and Q ~ r" 1 (Q ~ r _1//2 ) for a slow linear quench ending 
in the gapless (gapped) phase. We provide a qualitative explanation of such scaling. 

PACS numbers: 75.10.Jm, 05.70.Jk, 64.60. Ht 
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I. INTRODUCTION 



Non-equilibrium dynamics of quantum systems near 
quantum critical points has been a subject of intense 
study in recent years 1 ' 2 . During such dynamics, a quan- 
tum system passes from one gapped phase to another via 
time evolution of a Hamiltonian parameter A with a rate 
t -1 and an exponent a (A(t) = Ao|£/ r | a Sgn(i), where 
Sgn(x) = 1(— 1) for x > (<)0) through an intermediate 
quantum critical point at A = 0. At the critical point, the 
energy gap vanishes as A(fc) ~ \k\ z where z is the dynam- 
ical critical exponent. Thus the dynamics becomes non- 
adiabatic around a region near this point and the system 
fails to remain at the instantaneous ground state lead- 
ing to formation of defects 3-9 . The density of these de- 
fects (n) and the residual energy produced in the process 
(Q) scale with universal exponents: n ~ T -vda./(zua+i) 
and Q ~ T -{d+z)va/(zva+X) ^ wnere v j s the correlation 
length exponent and d is the system dimension 5,7 . It is 
well-known that scaling laws do not change if the dy- 
namics terminate at the critical point 8 . All of the above- 
mentioned studies apply to isotropic critical points where 
the scaling of the energy gap with the momentum is de- 
scribed by a single exponent z. Recently, the anisotropic 
Dirac model with an anisotropic critical point is studied 
and it was shown that one needs multiple exponents to 
describe the scaling of the energy gap 10 . However such 
studies have not been carried out in the context of the 
Kitaev model and generic expressions for the scaling laws 
for n and Q for such critical points in arbitrary dimen- 
sions have not been provided. Also, the effect of disorder 



on defect production in models, where the Harris crite- 
rion allows for the existence of a sharp quantum phase 
transition, has not been studied so far 11 . 

In this work, we study several aspects of non- 
equilibrium slow dynamics in the vicinity of both 
anisotropic critical points and critical points in the pres- 
ence of disorder with specific focus on the 2D Kitaev 
model which provides an explicit realization of both the 
cases. First, we derive a generic model-independent ex- 
pression for the scaling of n and Q for such dynamics 
which takes a d-dimcnsional system from a gapped phase 
to the vicinity of an anisotropic critical point. We con- 
sider a scenario where the energy gap vanishes as kf 
for m momentum components (i = l..m) and as kf for 
the rest d — m components (i = m + l..d) with z' > z at 
the critical point and show that the time-evolution of the 
Hamiltonian parameter A(t), which brings the system at 
the critical point at t = 0, leads to novel scaling laws for 
n and Q: 



— [m-\-(d—m)z / z']ya/ (zva-\-l) 
J ' ? 

- [{m+z) + {d—m)z/ z']i/a/ (zva+1) 



(1) 



Our results reproduce their well-known counterparts for 
the isotropic case (z = z') as special cases. We also show, 
by exact analytical solution for linear time evolution 
(a = 1), that the two-dimensional (2D) Kitaev model, in 
the absence of disorder, provides an explicit realization 
of the scaling laws mentioned above with d = z' = 2 and 
m = v = z = 1 leading ton~ r~ 3 / 4 and Q ~ t~ 5 / 4 . We 
also corroborate the scaling laws mentioned above by nu- 
merical studies of the Kitaev model for arbitrary power- 
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FIG. 1: Schematic representation of the Kitaev model on a 
honeycomb lattice. The bonds Ji, J2 and J3 shows nearest 
neighbor couplings between x, y and z components of the 
spins respectively, n represents the position vector of the 
midpoint of each vertical bond (unit cell). The vectors M\ 
and M2 are spanning vectors of the lattice. In the Fermionic 
representation of the model, the Majorana Fermions an and 
bn sit at the bottom and top sites respectively of the vertical 
bond with center coordinate n as shown. 



law time evolution. Second, we compute all independent 
multispin correlation function of the Kitaev model sub- 
sequent to a slow linear ramp which takes the system 
from a gapped phase to the vicinity of the anisotropic 
critical point, demonstrate their anisotropic nature, and 
discuss their spatial dependence. Third, we study non- 
equilibrium slow linear dynamics of the disordered Kitaev 
model where disorder is introduced via random choice 
of the fields Dn in the Fermionic representation of the 
model, and show, by explicit numerical calculation, that 
the defect production for such a dynamics obeys a differ- 
ent scaling law compared to its disorder free counterpart: 
n ~ r -1 / 2 and Q ~ r _1 (Q ~ r -1 / 2 ) for a quench ending 
in the gapless (gapped) phase. We provide a qualitative 
explanation for such defect production. 

The organization of the rest of the work is as follows. 
In Sec. II, we discuss scaling laws for defect density and 
residual energy for dynamics near an anisotropic critical 
point in the absence of disorder and show that the 2D 
Kitaev model constitutes an example of such a critical 
point. This is followed by Sec. Ill where we compute the 
equal-time correlation function of the 2D Kitaev model 
following such a dynamics and discuss its spatial struc- 
ture. In Sec. IV, we discuss defect production in the dis- 
ordered Kitaev model. Finally we provide a discussion of 
our results and conclude in Sec. V. 



II. ANISOTROPIC CRITICAL POINTS 

We begin with the study of slow dynamics in the 
Kitaev model 12 ' 13 . The Hamiltonian for this model, 
schematically represented in Fig. 1, is given by 

h k = ( j i t u t !+u + j 2 t !-iAi + j 3 t Ii t Ii+J> 

j-M— even 

(2) 

where fji — ( T jii T jiT T ji) denote Pauli matrices at the 
site of the honeycomb lattice, J±, J 2 , and J3 rep- 
resent nearest-neighbor couplings between x, y and z 
components of the spins respectively. It is well-known 
that Hk can be represented in terms of Fermionic 
fields by a straightforward Majorana transformation: 

«ii = (rii=-oo t !i) r ji for even 3 + 1 and h = 

(nt-00 T i) T ji for odd 3 + l13 - This leads t0 the 
Fermionic Hamiltonian 

Hf = i t Jl b ^ a H-m + 1/2 6 « a «+A? 2 

ft 

+J 3 Dn bjtan], (3) 

where n = m + ( i + %fj ri2 denote the midpoints 
of the vertical bonds. Here ni, 712 run over all integers so 
that the vectors n form a triangular lattice whose vertices 
lie at the centers of the vertical bonds of the underlying 
honeycomb lattice. The Majorana Fermions an and bft 
sit at the bottom and top sites respectively of the bond 
labeled n. The vectors M x = ^-§j and M 2 = ^+§J 
are spanning vectors for the lattice, and Dft can take the 
values ±1 independently for each n. The crucial point 
that makes the solution of Kitaev model feasible is that 
Dn commutes with Hf, so that all the eigenstates of Hf 
can be labeled by specific values of D^. It is well-known 
that the ground state of the model corresponds to = 1 
on all links 12 . 

For Dn = 1, Eq. (3) can be diagonalized as 

k 

where ipt — (at, bt) are Fourier transforms of an and 

bn, the sum over k extends over half the Brillouin zone 
(BZ) of the triangular lattice formed by the vectors n, 
and Hj: can be expressed in terms of the Pauli matrices 
a 1 in particle-hole space as 

H % = 2[Jisin(fc-Mi)- J 2 sin(fc-M 2 )]cr 1 

+2[J 3 + Ji cos(fc • Mi) + J 2 cos(fc • M 2 )]er 2 . (5) 

The spectrum consists of two bands with energies £? = 
±E^ 13 , where 

Ej: = 2[{Jisin(£- Mi) - J 2 sin(fc • M 2 )} 2 

+{ J 3 + Ji cos(£ • Mi) + J 2 cos(fc ■ M 2 )} 2 ] 1 / 2 (6) 
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For | Ji — J2I < J3 < J1 + J2, the bands touch each other, 
and the energy gap A; = Ei — EZ vanishes for special 

K k k 

values of k leading to a gapless phase. In particular we 
note that for J\ = J 2 = 1 and J3 = 2, the gap vanishes 
at k c — (27t/-\/3, 0) and around this point Aj ~ k y and 
A ^ <~ k 2 . Thus this critical point constitutes an example 



the limits of integration can now be safely extended to 
infinity. To compute this integral, we note that around 
k = k c , a s - k ~ 3Sk y and g g - k — 3((5fc 2 + 3<5fc 2 )/4. Thus 
a redefinition of variables 8k x — > 5k' x — Sk^ 1 ^ 4 and 
Sky — » Sky = Sk y T X l 2 allows us to extract the r depen- 
dence of the defect density 



of an anisotropic critical point with z = m = 1 and d = „ „ 

2. We note that such an anisotropic scaling occurs n ~ / d5k x dSk v p s ~ k , <~ r~ 3//4 / dSk' x dSk' y p s p. (10) 

nv nnn-7Prn vfllnp nf .7^ and at = ( -4- ^ ^ 



1 and <i = 



for any non-zero value of J\ and J 2 at J 3 = (J x + J 2 ). 

We now consider a dynamics in this model Js(t) = 
(Ji + J2 — Jt/r) from t = —00 to t = at a fixed rate 
1 /r which brings the system from a gapped phase to the 
anisotropic critical point at k c . Although this quench 
problem can be solved for any Ji and J 2 , we shall fix 
Ji = J 2 = J for simplicity and scale all energies (times) 
by J (h/J) in the subsequent analysis. This choice does 
not change the scaling properties which we seek. Also, 
to study the time evolution of the system, we note that 
after an unitary transformation U — exp(— zcr 1 7r/4), we 

obtain H F = EcV'r #pV>n where HL = UHrW is given 

— K k k k *• K 

by 



H', = 2[(g % - t/r)a 3 



(7) 



where 



sin(fc • Mi) — s\n(k ■ M 2 ) and g^ = 2 + cos(A: • 

Mi) + cos(fc • M 2 ). Hence the off-diagonal elements of 
H'~ remain time independent, and the quench problem 

reduces to a Landau- Zener problem for each k. 

The state of the system after the quench at t = can 
be found by solving the Landau- Zener problem at each 
k with the initial condition ip^{t = —00) = |1) = (0, 1) T 

for all k. After some algebra, one obtains for a given k 
and at t = 14 



e -7ra|r/4/ e 3 l7 r/4 



+a^D^i(^)\0)), 



(8) 



where 



2ig^y/r exp(— iir/4), /i£ = —lair and are 
parabolic cylinder functions. The excited state at t = 0, 
solved by diagonalizing HL(t — 0), yields, for a given 

= u E i - 2g ^ + 2a k\ Q ))/ v v where v k = 

[(Ei - 2g~) 2 + 4a|] 1 /2. Thus the probability of defect 
formation, given by p^ — \(^i\ip k ) d \ 2 , can be obtained 



as 



n = 



4ale 
k 



-7rair/2 



Ei 

k 



(9) 



Since r is large for slow dynamics, the contribution to 
the defect formation comes from a small region near the 
critical point where Ag is sufficiently small for k ~ k c . 
The density of defects can be thus estimated by expand- 
ing pj: about k = k c : n ~ / dSk x dSk y P% = % c+ f k , where 



A similar analysis can be carried out for computa- 
tion of residual energy Q = (2ir)~ 2 J d 2 kp^A^. Here 

we note that near the critical point k — k c , Ag ~ 
\hSkl + 9(6kl + 3(5fc2) 2 /16 and thus scale as t" 1 / 2 . 



Thus one obtains 



Q — I dSk x dSky Af k Pf k , ~ t 



-5/4 



dSk' x dSk' y A spPsP . 



(11) 



Eqs. (10) and (11) show that n ~ r~ 3 / 4 and Q - t~ 5 / 4 
at the critical point. These scaling laws do not conform to 
the predictions of earlier works on defect production dur- 
ing passage through isotropic quantum critical points 5 or 
critical surfaces 13 ; their origin lies in the anisotropic scal- 
ing of Sk x and Sk y with the quench time r. 

To generalize these results for arbitrary d-dimcnsional 
anisotropic critical points, where the energy gap A^ <~ kf 
for m directions and <~ fcf for d~m directions, we provide 
a simple phase space argument as first proposed in Ref. 
4. We consider a general power-law quench with A(t) = 
Ao|i/r| Q Sgn(t) which starts at t = —00 and reaches the 
critical point at t = 0. We first note that the adiabaticity 
condition breaks down when the rate of change of the 
energy gap become equivalent to the square of the gap: 
dA^/dt > A|. Since A £ - X zl/ \t/r\ zl/a , we find that the 
time spent by the system in the non-adiabatic regime is 



zvcuj (zva 



+1 K The scaling of the energy gap 



-ZVOij (ZVCL+1) 



given by t ~ r 

in this regime can thus be written as A^ 
The phase space for defect production is given by f2„ 
k\..kd- Since A^ <~ kf for i = l..m and fcf for i 
m + l..d, we finally obtain 



n * 



-{m+{d—m)z/z')va/ (21/a+l) 



(12) 



A similar argument can also be presented for the residual 
energy. We note that for z < z', the leading behavior of 
the energy gap near the quantum critical point, where 
the defects are produced, is A^ <~ kf for 1 < i < m. 
Thus the phase space for the residual energy production 
is Oq ~ A^ki-.kd leading to a scaling of Q as 



Q 



-\{m-\-z)-\-{d—m)z I z']vct / (zva+1) 



(13) 



We note that the scaling laws, Eqs. (12) and (13), repro- 
duce their isotropic counterparts for z — z' leading to 

n ^ T -d Va /(zva+l) and Q _ T -(d+z) Va /(z Va +l)5,7,8_ AlsQ; 

the scaling of the Kitaev model for linear time evolution 
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elaborated in this work is reproduced for d = z' = 2, 
and z = v = a = 1 leading n ~ r -3 / 4 and Q ~ r -5 / 4 . 
Moreover, we note that the scaling of defect density for 
a linear quench through a gapless surface can also be ob- 
tained from Eq. (12) by noting that for such quenches the 
energy gap depends only on the m momenta components 
orthogonal to the d—m dimensional gapless surface. This 

1 1 z 1 

can be represented by putting z' — >• oo (since k\\ ~ A- ) 
leading to the scaling law n ~ T -mu/{zu+i)i3_ Thug Eqg _ 
(12) and (13) reproduce all earlier results on defect pro- 
duction for slow dynamics across quantum critical lines 
and surfaces as special cases. Finally, we would like to 
point out that the maximum values of these exponents is 
2 which can be obtained by similar considerations as in 
the cases of isotropic critical points 8 . 

To verify these scaling laws, we now study non-linear 
power-law dynamics in the Kitaev model numerically. To 
this end, we again restrict ourselves to J\ = J% = 1 and 
evolve J 3 (t) = (2 - |t/r| a Sgn(t)) for -oo < t < so that 
the anisotropic critical point is reached at t = 0. The 
corresponding time-dependent Hamiltonian is given by 
H(k; t) = J2k 4 [Gfc ~ Wr\ a Sgn(t))a 3 + a^a 1 ] ^. We 
solve the time-dependent Schrodinger equation idtip^ = 
H(k;t)ip(k) numerically for each k, compute pt, and 
use it to obtain the defect density n = J d 2 kp^ and 
Q = f d kAjzp? numerically as a function of r and a. 
The plots of n and Q vs t are shown in Fig. 2 for sev- 
eral representative values of a. The lines in the figure 
indicate the power laws expected from Eqs. 12 and 13 
(n ~ r- 3a /[ 2 ( Q+1 )l and Q ~ r -6«/P(«+i)l) for d = z' = 2 
and z = v = m = 1. The agreement between the nu- 
merical and theoretical results corroborates the scaling 
theory proposed in this work. 



III. CORRELATION FUNCTION 
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FIG. 2: Numerical results on the defect density n and residual 
energy Q. The time-dependent Schrodinger equation is solved 
in the momentum space for systems with size up to 512 x 
512 unit cells. The parameter a specifying the evolution of 
J3 = 2 — |t/rj Q Sgn(t) is chosen as a = 1, 3 and 5. The lines 
indicate the power laws expected from Eqs. (12) and (13) for 
d=z' = 2smdz = v = m=l, n~ r - 3a/[2(a!+1)] and 
Q ~ T - 5a /l 2 ( a + 1 )\ T ne agreement between curves obtained 
numerically and the corresponding power laws is remarkable. 
In all plots, t varies from an initial value ii n = — 3r to a final 
value t{ — 0. 



In this section, we compute the independent correla- 
tion function for the Kitaev model for linear time evo- 
lution. Since the model can be represented by free 
Fermions, it is easy to see that the only non-zero in- 
dependent correlators are those between free Fermions 
which are given by 

(Of) = i{bftOri+i>) 

An 

= ^E[(&k)exp(zfc.r)-h.c], (14) 



k 



where (..) denotes expectation value with respect to a 
direct product of states involving k only, h.c. denotes 
hermitian conjugate, and N s is the number of sites. After 
an unitary transformation U = exp(— ia tt/4), we find 

2 

N 



pointed out in Ref. 13. For r = 0, (Op) represents corre- 
lations between z components between spins at the end 
of the vertical bond whose midpoint is denoted by n. For 
r 7^ 0, it represents correlation between product of multi- 
ple spin operators which begins with t x or t v on a b or a 
site at ft — (j, I) and ends with t x or t v on an a or b site 
at n + r = (j', V) with a string of r 2 operators living on 
sites in between. Note that the Fermionic representation 
in terms of free Fermions with _D,-j = 1 ensures that these 
multispin correlation functions are the only non-zero in- 
dependent spin correlation functions of the model. 
The ground state with a fixed k for J3 = 2 is given 



by \ip%. 



- 2 9i )\l) + 2a^\0))/VZ, where V. = 
21 1 / 2 . Noting that |0) and |1) are basis 



(Of) = -— ^(^- t [-cos(A:-r)fT 3 +sm(A:-r)(T 1 ]'i/;^). (15) of the ground state as 



of ipL, one finds, from Eq. (15), the correlation function 



The interpretation of these correlation functions in terms (Op) 
of the original spin degrees of freedom have already been 
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FIG. 3: Plot of the correlation function {Of) de ^ ect as a func- 
tion of m and «2 for linear quench of J3 / J\ from 5 to 2 with 
Ji = J2 = 1 and r = 5. 



sin(fc • r) 



(16) 



where A = 47T/3V3 is the area of half of the Brillouin 
zone. As for the state after quench, a straightforward 
calculation using Eq. (8) shows 

(O r -) d = Sn + j [ d 2 ke-™l T ' 2 a % ^[2a^ 



xD^i/g) + c.c.j sin(fc • r) 



(17) 



Note that (Of) d reduces to (0?) G with r -> 00. The 
correlation between defects induced by non-adiabatic 
quench dynamics can be captured by the deviation of 
{Op) d from (Op) G . Thus we define the defect correlation 
function by 

(O r -) dofcct = (Of) d - (0 P ) G . (18) 

The nature of the spatial dependence of the defect corre- 
lation function for slow dynamics (large r) can be qual- 
itatively understood for J\ — J 2 from Eq. (17). To this 
end, we first separate the contribution to (Op) d which 
comes from around aj; ~ 1/t from those coming from 

other regions in the fc-space. For estimating the latter 
contribution, we consider t > 1 so that for a^g^ 7^ 0, 

l^fel' \ v %\ ^ 00 ■ We * nen n0 * e that * ne following identi- 
ties for D holds in the limit b — > 00 with arbitrary ratio 
a/b 



-■Kb 2 / 4 



D. 



e -T&74£)_. 62 ( ae W4) „ cos (6)e- ir >, 



-i(rj+ir/4) 



(19) 



where 8 and 77 are defined through the relations 



cos(0)(sm(0)) = V /[l + (-)a/(2 A /6 2 + a 2 /4)]/2, 
7? = -6 2 /2 + b 2 ln(a/2 + V^ 2 + a 2 /4) 

+a v / 6 2 + a 2 /4/2. (20) 



FIG. 4: Plot of the peak positions of {Or) defect in the m - n 2 
plane for J\ = 1 and several representative values of J2. For 
each of these cases, the quench starts at J3/J1 =5 and ends 
at the anisotropic critical point. Note that the axis of 112 is 
upside down. 



Identifying b = ol^^/t and a = 2g^r and substitut- 
ing Eqs. (19) and (20) in Eq. (17), we find, after some 
straightforward algebra, that the integrand of Eq. (17) 
reduces to that of Eq. (16) for all k except those for 
which cxts/r — 1. Thus, one finds that in this limit, the 
main contribution to (CV) dcfcct comes from around the 
line ag ~ 1/y/r. For large r, this is infinitesimally close 
to the line sin(fc • Mi) = sin(fe • M2). In this region of k 
space, |D Ms _ 1 (!/g)| 2 ~ |A-i(^)| 2 which, for large r and 
J3 = 2, is a sharply peaked function for c/^ ~ which oc- 
curs at k = k c . Also, for k ~ fc c , it can be easily checked 
that [D*_ 1 (^)D^(^)+h.c] « iD^.r^)! 2 , so that 

the major contribution to (CV) dofoct comes from the co- 
efficient of the cos(A: • r) in the integrand. Using these 
observations and expressing r = + n.2/2), 3^2/2), 

one can estimate the spatial dependence of the correla- 
tion function in the same line as in Ref. 13. In particular, 
the maxima of the correlation function is expected to oc- 
cur along the maxima of cos(fc c ■ r) i.e. along the line 
ni + 712/2 = in the ri\ — n-i plane. Away from this line, 
as shown in Ref. 13, (CV) dofect is expected to decay ex- 
ponentially as a function of r with a characteristic decay 
length ~ y/r. 

A plot of (Op) dcicct as a function of n\ and 712, obtained 
by numerical evaluation of Eq. (18) are shown in Fig. 
3 for J\ = J2 = 1, corroborates the above-mentioned 
discussion. We find that (0?} deiect peaks along the n\ — 
—712/2 line and decays to zero as we move away from this 
line. The decay length in the n\ — n<i depends on r; for 
larger r we have a sharper decay. The slope of the line 
along which (Op) dclcct peaks in the n\ — 122 plane changes 
with J1/J2 since k c depends on this ratio. This can be 
seen from Fig. 4 which plots the position of the peaks of 
the correlation functions for several representative values 
of J\j ' J%. The analysis of the preceding paragraph can 
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be easily extended to these cases in the same line as in 
Ref. 13 and is found to match the numerical results for 
all Ji/J 2 . 



IV. DISORDERED KITAEV MODEL 

In this section, we study the dynamics of Kitaev model 
given by Eq. (4) with a random configuration of D^, 
namely for random assignment of values ±1 to the link 
variables D^. The dynamics is incorporated in the form 
of a power-law evolution of J3 as in Sec. II. 

The Hamiltonian (Eq. (3)) can be expressed using the 
real-space Fermion operators an = ^(bff—iaft) at position 
ft as 

h f = J2 Ji ( a « + 4) ("s-fA - a l 

n 

+J 2 (a H + at) (*n+M 



M 2 n+M 2 , 

J 3 D H (1 - 2ala a ). (21) 



The first two terms represent hopping and pair-creation 
and annihilation of the Fermions while the third term in- 
duces a random local potential. For J3 J\ %, the third 
term dominates and the ground state of the system is 
composed of localized states of the Fermion. In contrast 
for J3 = 0, the ground state is clearly delocalized. We 
now show numerically that a quantum phase transition 
takes place in between these two limits at J3 = J^, c - Note 
that the existence of a sharp transition in the presence of 
the disorder is consistent with the Harris criteria vd > 2 
since for the Kitaev model d = 2 and v = 1. 

The Hamiltonian, Eq. (21), is written in a 
quadratic form as H = ip^Mtp with ifi = 

' a ff 2 ' ■ ■ ■ ' a «™ ' a "i ' ' • -i^jv)) where N is the 
number of vertical bonds (unit cells) in the system, and 
M is a 2N x 2N matrix given by 



M 



with 



-Mi 



A B 
B T -A 



(22) 



fii-\-M2 : ni fti ,rti~{-M2 



Hi— Mi ,fii 



2J 3 D n 

-B 



Ui ,fli —Ml 



-Ji, 
-Jl, 



(23) 



All other elements of A and B are zero. The matrix M 
is diagonalized by a unitary matrix, 



U 



u v 

* 

V u 



(24) 



as WMU = D, where D is a diagonal matrix. We note 
that the form of M necessitates that if e M is an eigenvalue 
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FIG. 5: Probability distribution of excitation gaps at J3 = 
1.5. 10000 instances of {-Dh} are generated and for each of 
them we obtained the excitation gap 2ei by numerically di- 
agonalizing the matrix M. 



of M, so is — e M . We hereafter suppose e p > and choose 
U so that e M (/x = 1, 2, • • • , N) enter upper half diagonal 
elements of D. Defining a fermion operator as 



7p =5Z u «>M a « + v k» a fi> 



the diagonalized Hamiltonian is written as 



(25) 



N 



(27^ - !)■ 



(26) 



The ground-state energy is given by E g 
The energy gap from the ground state to the first excited 
state is thus given by A = 2ei where t\ is the smallest 
positive eigenvalue. 

With this observation, we now compute the gap A nu- 
merically for finite sizes and obtain the distribution of 
gaps by changing the configuration of \Dn\- We find 
that the property of the distribution of gaps is quali- 
tatively different for J3 < 1.5 and J3 > 1.5. Let us 
first consider the case with J3 = 1.5 for which the dis- 
tribution of the gaps is shown in Fig. 5 for several sys- 
tem sizes. The shown distribution allows a Gaussian fit: 
iV(A) = (2Tra 2 y 1/2 e- < - A - A ^ 2 ^ using the average A 

and the variance a 2 = A 2 — A , where ~ stands for the 
average over the random configuration of {D^}. Figure 6 
shows the size scaling of A for several J3 > 1.5. We find 
A scales linearly with 1/L. Since the variance of gaps 
tends to vanish for L — ¥ 00, one can estimate the gap in 
the thermodynamic limit by extrapolating the fit- 
ting line of A for 1/L — > 0. Such a behavior of Aoo is 
to be contrasted with that for J3 = 1 as shown in Fig. 
7. For J3 = 1, we find that the probability distribution 
of gaps scales as A ~ 1/L 2 as seen from the collapse of 
the data for several system sizes (Fig. 7) . The difference 
in behavior of A can be further understood by plotting 
Aoo for several values of J3. This is shown in Fig. 8. We 
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FIG. 6: Finite size scaling of the average of gaps for J3 = 1.5, 
1.6, 1.7, 1.8, 1.9, and 2.0. We find that the average of gaps is 
scaled by l/L, where L is the length of the system (iV = L 2 ). 

find that for J 3 <, 1.4 the gap vanishes, while it increases 
linearly with J3 for J3 > 1.5. The position of the crit- 
ical point J3 iC can therefore be estimated to be around 
1.5. Moreover, the gap A increases as Aa | J3 — J3 jC | for 
J3 > 1.5 leading to zv = 1 for the transition. 

Having established the presence of a quantum critical 
point in the disordered Kitaev model, we now study the 
dynamical behavior during slow non-adiabatic linear time 
evolution J^(t) — —Jt/r which takes the system from a 
gapped region ( J 3 = 5) either to a gapless region ( J 3 = 0) 
or to a gapped region passing through the gapless region 
(J3 = —5). In order to obtain quantities of interest, 
we switch to the Heisenberg picture 15,16 and introduce 
the time-evolution operator U(t), = U (t)\^ (t in )) , 

where t- m denotes the initial time. The operator in the 
Heisenberg picture is denoted by = w (t)ctfiU{t). 

Computing the commutator of cm and Hp expressed by 
Eq. (21), the Heisenberg equation of motion for a^i is 
given by 

= E (Mn°%(t) + B^afitj) . (27) 

n 

We define matrices u^ lV (t) and v^^it) by an expan- 
sion of a^i (t) by 7„.in, operators which diagonalize the 
Hamiltonian at initial time t m (see Eq. (26)): 

"*(*) = E (<wW7,,m + «&,„(t)7ita) ■ (28) 

V 

Substituting this expansion for a H 's in Eq. (27), one ob- 
tains equations of motion for u^ v (t) and v*^ v (t)\ 

n 
n 

The initial conditions for u^^if) and u^ i( ,(t) are written 
as Mm^(tin) = Um,i/,in and ^^(fin) = WrS^.in: where u ia 



FIG. 7: The probability distribution of gaps at J3 = 1.0. 
Horizontal axis is the excitation gap multiplied by L 2 . The 
curves with different size almost collapse, meaning that the 
distribution of gaps is given by a function of AL 2 and the gap 
vanishes as l/L 2 with increasing L. 

and Vi n are block matrices of U diagonalizing M at initial 
time. To obtain the expressions of n and Q at final time 
tf, we introduce notations itf, Wf, e Mj f, and 7 M .f so that 

H F {U) = E M e M,f( 2 7if7M,f - !)> where 

7^f = X! u ^,f a ™ + v a,/*,f a ff ( 31 ) 

n 

The density of excitation n and the residual energy Q 
can now be defined by 

1 W 

Q = (*(tt)\H F (tt)\*(tf)) - E g . (32) 

Next, we switch to the Fermion operators a from 7f using 
Eq. (31) and shift to the Heisenberg representation. Sub- 
stituting the expansion Eq. (28) in Eq. (32), one obtains 

n = Itr [(vfuitf) + ujv{t()) («t(t f )uf + ^(tfK)] , 

iV 

Q = 2^e M , f 

x [(« f T u(tf) + u^(t f )) (« f (*fK + «t(tfK)] Mi|i . 

(33) 

First, we present numerical results for two cases of the 
quench. For each value of the quench time r, simula- 
tions were carried for 16 different configurations of {D r i} 
for obtaining a large enough sample set for disorder av- 
eraging. Figure 9 shows the disorder averaged values of 
density of excitations and residual energy as a function 
of r after the time evolution. The results of simulation 
suggest that for large r, the density of excitation n and 
residual energy Q scale with r as 

n ~ r- 1 / 2 , Q~t~\ (34) 
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FIG. 9: Scalings of the density of excitations and residual 
energy after a quench of J3 from 5 to —5 and from 5 to 0. 
The density of excitations is scaled as n ex ~ 7- -1 / 2 in both 
cases. The scaling of residual energy is Q ~ r _1 when J3 
stops at 5 and Q ~ t~ 1//2 when J3 stops at 0. Simulations 
are carried out for systems with 16 x 16 unit cells. The average 
is taken over 16 configurations of {Dh}. 



for an evolution ending inside a gapless phase and 

n ~ r- x /\ Q ~ r" 1 / 2 , (35) 

for that ending in a gapped phase after passing through 
the gapless phase. We note that these scaling laws are 
different from those obtained for uniform Dft . 

A qualitative explanation of such scaling laws for n 
and Q can be obtained as follows. We recall that for 
dynamics in critical systems without disorder, the con- 
dition for diabaticity is given by ^ > A 2 (Ref. 1). In 
generic second order quantum phase transition with crit- 
ical exponents z and v, one can write A ~ \ zv where 
A is quenched with a rate 1/r. This yields standard 



FIG. 10: Density of states of quasi-particles at Ji = J2 = 
J3 = 1 in gapless phase. The quasi-particle energies e M are 
computed for systems with 32 x 32 unit cells and the histogram 
of them is obtained. The bin is set at 0.1. The average is 
takes over 10000 instances of {-Dh}- The density of states 
with small energy takes a finite value and fluctuates. The 
amplitude of fluctuation is comparable with the error bars. 
This result suggests D(e) is a constant when e is small. 



expressions 



-zv j (zv-\-l) 



From this, one can es- 



timate the scaling form of the density of excitations and 
the residual energies to be 



D(e)de, Q 



(A f + s)D{e)ds, (36) 



where D(e) is the density of states of quasi-particles near 
the critical point or gapless region and Af is the final 
excitation gap when the quench stops. Note that Af = 
for a quench ending in the gapless region. Typically, the 
density of states at the critical point or in a gapless region 
is given by D(e) ~ e p for some non- negative exponent 
p. Using this, one may obtain scaling of the density of 
excitation and the residual energies as 



^p+l ^ T -{p+l)zu/(zv+l) 



-(p+2)zv/(zv+l) 



Q ^ A p+2 - T 



(37) 



where in the second line we have assumed that Af = 0. 
For finite Af,n and Q scales according to the same power 
law. 

To obtain the scaling of the gap, we need to obtain the 
value of p. To this end, we plot the density of states for 
a finite-sized system with 32 x 32 unit cells in Fig. 10. 
The plot suggests that the density of states is a constant 
at least at low energies e/J< 0.5 implying p = for the 
critical modes. We have checked that this holds for other 
system sizes as well. Moreover, numerical studies shown 
in Fig. 8 leads to zv = 1. Using these facts, one obtains 



'1/2 



.-1/2 



gapless phase 
gapped phase. 



(38) 



The scaling laws in Eq. (38) can be also obtained by 
another argument. To elucidate this, we show, in Fig. 
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Note that these arguments do not depend on whether 
the quench ends inside the gapless phase or not since for 
slow dynamics the defects are produced mostly during 
the passage through the gapless regime. Using the fact 
that p — for these systems, a similar analysis for Q 
reproduces the results of Eq. (38). 



DISCUSSION 



FIG. 11: Low- lying positive eigenvalues of M as a function 
of J3 with a fixed configuration of {Da}- The system is com- 
posed of 8 x 8 unit cells. 



11, the low- lying energy spectra of quasi-particles of the 
model as a function of J3 for finite sized system (8x8 
unit cells) for a single configuration of {-D,i}. Since there 
is a finite gap for all values of J3, an adiabatic evolu- 
tion do not lead to quasi-particle excitation. For non- 
adiabatic processes, the most probable excitation occurs 
around the avoided level crossing with minimum energy 
gap shown with a blue arrow in Fig. 11. We denote 
the corresponding energy gap by A;. The probability 
of excitation is well approximated by the Landau-Zener 
formula: e~ cAlT , where c is a constant factor determined 
by the slope of the excitation gap around A/. Next, we 
recall that the distribution of excitation gaps A for fixed 
J3 inside the gapless phase scales as 1/L 2 . Hence the 
distribution of A; is also a function of A/L 2 . Thus the 
probability distribution function of A; can be written 
as P{u) with u — A;i 2 . Assuming that the factor c 
is independent of L and A;, the averaged probability of 
excitation n is given by 



duP(u)e~ cA ' 
H(r/L 4 ). 



duP(u)e- cu2T/Li 



(39) 



From this, one can obtain a length L £ that yields av- 
eraged probability of excitation II = e for a given r: 
1/4 



(<0 



Nr. 



For sufficiently small e, 

regarded as the average size within which a single excita- 
tion is expected to occur. The density of these excitations 
is thus estimated by N s as 



1 

nZ 



n-i(e) 



1/2 



-1/2 



(40) 



In conclusion, we have shown that the Kitaev model 
constitutes an example of a two-dimensional model with 
an anisotropic critical point. We have also demonstrated 
that the presence of such an anisotropic critical point 
leads to novel scaling laws defect density and residual 
energy during slow power-law dynamics which takes the 
system from a gapped phase to the vicinity of such a crit- 
ical point. We have generalized our results for such scal- 
ing laws for d-dimensional systems with such anisotropic 
critical point. Further, we have computed all indepen- 
dent correlation functions of the Kitaev model in the 
Fermionic representation after a slow linear ramp which 
brings the system to the vicinity of an anisotropic criti- 
cal point. We have charted out the spatial dependence of 
the correlation function and discussed its relation with 
several multiple spin correlators of the model. Finally, 
we have studied the non-equilibrium slow dynamics of 
the disordered Kitaev model where disorder is introduced 
via random configuration of Dfi in its Fermionic repre- 
sentation. We have shown numerically that the defect 
density n, generated during a slow linear ramp from a 
gapped phase of the model to either a gapless phase or 
to another gapped phase through a gapless region, scales 
as r -1 / 2 . In contrast, the residual energy Q scales as 
t -1 / 2 (t -1 ) for similar dynamics ending on the gapless 
surface (gapped phase after passing through the gapless 
surface) . We provide a qualitative understanding of such 
scaling laws to back up our numerical results. We note 
that there has been suggestions of experimental realiza- 
tion of the Kitaev model using ultracold atomic system 17 . 
In the event of such a realization, the simplest experi- 
mental test of our theory would involve measurement of 
defect density n following a slow ramp. Such experiments 
has recently been performed for standard ultracold boson 
systems 
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